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Abstract 



The symmetric multistep methods developed by Quinlan and Tremaine (1990) are shown 
to suffer from resonances and instabilities at special stepsizes when used to integrate nonlinear 
equations. This property of symmetric multistep methods was missed in previous studies that 
considered only the linear stability of the methods. The resonances and instabilities are worse for 
high-order methods than for low-order methods, and the number of bad stepsizes increases with 
the number frequencies present in the solution. Symmetric methods are still recommended for 
some problems, including long-term integrations of planetary orbits, but the high-order methods 
must be used with caution. 

1 Introduction 

In many physical problems one has to solve second-order differential equations of the type 



where x{t) is the position at time t and f(x, t) is the force, assumed to be independent of the velocity. 
Such an equation or set of equations can be solved efficiently by a multistep method 



where x n is the computed position at time step n and h is the stepsize. A popular class of such 
methods is the the Stormer-Cowell class, for which a k = 1, a k -\ = ~~ 2, ctk—2 = 1> and a^s = 
■ ■ ■ = ctQ = 0, with the Stormer method being explicit {[3 k = 0) and the Cowell method implicit 
(Pk 0)- Stormer-Cowell methods have often been used for long-term integrations of planetary 
orbits (see Quinlan and Tremaine 1990 and references therein). But the Stormer-Cowell methods 
suffer from a defect, sometimes called an orbital instability, when the stepnumber k exceeds 2: if a 
Stormer-Cowell method with k > 2 is used to integrate a circular orbit, the radius does not remain 
constant, and the orbit spirals either inwards or outwards (the direction depends on k). This defect 
was recognized by Gautschi (1961) and Stiefel and Bettis (1969), who proposed modified multistep 
methods for orbital integrations. Their methods require a priori knowledge of the frequency of the 
solution, however, which is usually unknown or, at best, known only approximately. 
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Lambert and Watson (1976) showed that the orbital instability of the Stormer-Cowell methods 
can be avoided by choosing the coefficients of a multistep method to be symmetric, so that 

at = a k -i, Pi = f3 k _i, i = 0,...,k. (3) 

Lambert and Watson analysed in detail the application of symmetric methods to the linear test 
equation 

x"{t) = -cA(i), (4) 

and showed that if ui 2 h 2 lies within an interval (0, Hq ), which they called the interval of periodicity, 
the solution is guaranteed to be periodic. Quinlan and Tremaine (1990) extended the work of Lam- 
bert and Watson (1976) to derive high-order explicit symmetric methods suitable for the integration 
of planetary orbits, and compared these with high-order Stormer methods. The symmetric methods 
gave energy errors that did not grow with time, and position errors that grew only linearly with 
time, whereas the Stormer methods gave energy errors that grew linearly with time, and position 
errors that grew as the time squared. 

Soon after Quinlan and Tremaine's methods were published, however, Alar Toomre discovered 
a disturbing feature of the methods, an example of which is shown in Figure H. Panel (a) shows 
the maximum error in the energy of a circular Kepler orbit integrated with the 8th-order symmetric 
method of Quinlan and Tremaine and with the 8th-order Stormer method, plotted versus the stepsize 
used in the integration. The energy error decreases with the stepsize as ~ /i 9 , as expected for an 
8th-order method, and at most stepsizes the error from the symmetric method is much smaller than 
the error from the Stormer method. But there is a startling spike in the energy error from the 
symmetric method near a stepsize corresponding to 60 steps per orbit, a stepsize that is well within 
this method's interval of periodicity. Figure [2] shows the time development of the instability for a 
circular orbit integrated with 60 steps per orbit. The energy error grows exponentially for about 
the first 400 periods until it reaches a maximum valueQ of about 0.25. It decreases over the next 
few hundred periods to a minimum value of about 10~ 7 , and then oscillates between these extremes 
with a period of roughly 550 orbital periods. The longitude error (not shown in the figure) grows 
exponentially until it reaches a value of order unity and then stays at that level. 

The problems are worse for eccentric orbits, as shown in panel (b) of Figure |l| for an orbit with 
e = 0.2. The symmetric method is still much better than the Stormer method at most stepsizes, 
but now there are spikes in the energy error at a number of stepsizes. The spikes that appear at the 
stepnumbers (the number of steps per orbit) 90, 120, 150, etc. are similar to the spike at 60 steps 
per orbit: they are instabilities at which the error grows exponentially with time. The smaller spikes 
at stepnumbers that are multiples of 5 or 6 (54, 55, 65, 66, 70, 72, etc.) are resonances, and not 
instabilities, since at these stepsizes the energy error grows linearly with time, as shown in panel (b) 
of Figure 

The resonances and instabilities do not occur for the linear test equation (|]) that has been used 
in previous discussions of symmetric multistep methods. The present paper is written to explain 
their origin, to see if anything can be done to reduce their severity, and to decide if they render the 
methods unsuitable for problems like the long-term planetary integrations considered by Quinlan 
and Tremaine (1990). 

lr The value of the energy error at the instability depends on the formula used to compute the velocities. If the 
multistep equation is written in summed form (Henrici 1962; Quinlan 1994) and the velocities are computed from 
the summed accelerations, the maximum error is about 10 times smaller than shown in Figure ^, but apart from this 
difference the instability remains the same. In this paper the summed form was used only for the planetary integrations 
described in Section n. 
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Figure 1: Maximum fractional energy error during integrations of a Kepler orbit for 25000 periods using 
the method SY8 (solid lines) and the 8th-order Stormcr method (dashed lines), plotted as a function of the 
number of steps per orbit: (a) circular orbit; (b) eccentric orbit (e = 0.2). 

2 Symmetric multistep methods 

We start by reviewing symmetric multistep methods and some properties of them that are needed 
for explaining the resonances and instabilities. These properties are described in more detail by 
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Figure 2: Time dependence of the fractional energy error with the 8th-order symmetric method SY8. Panels 
(a) and (b) plot the maximum error observed in the previous period during integrations of Kepler orbits with 
(a) a circular orbit using 60 steps per period and (b) an eccentric orbit (e = 0.2) using 54.03 steps per period. 

Henrici (1962) and Lambert and Watson (1976). A survey of multistep methods and other methods 
for integrating second-order differential equations is given by Coleman (1993). Practical techniques 
for reducing roundoff error in long multistep integrations are described by Quinlan (1994). 

We consider fc-step multistep methods of the type (||), where without loss of generality we assume 
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otk = 1 and |ckq| + |/3q| > 0. With the multistep method we associate the linear difference operator 

k 

L [x(t);h] = [ujx(t + j-h)- h 2 fyx" '(t + j ■ h) . (5) 

3=0 

If x(t) has continuous derivatives of sufficiently high order then 

L [x(t);h] = C x(t) + dx'(t)h + ■■■ + C q x {q) {t)h q + ■■■, (6) 

where 

C q = ^(0 q a + ■■■ + k q a k ) - —L-(W~ 2 (3 + • • • + k q ~ 2 f3 k ) (7) 

(the second term on the right is absent if q < 2). The order p is the integer for which Co = 
••• = C p+ i = 0, C p+ 2 7^ 0. A method is said to be consistent if its order is at least 1, i.e., if 
Co = C\ = Ci = 0. 

To discuss the stability of multistep methods we introduce the polynomials 

p(z) = a k z k + ak-iz^ 1 H h«o, (8) 

a(z) = P k z k + f3k~iz k ~ l + ■ ■ ■ + Po. (9) 

A method is said to be zero-stable if no root of p{z) has modulus greater than one and if every root 
of modulus one has multiplicity not greater than two. A consistent method has p(l) = p'{T) = 0, 
so for zero-stability p(z) must have a double root at z = 1. This is called the principal root; the 
other roots are called spurious. A method is convergent if and only if it is consistent and zero-stable 
(convergence means essentially that x n — * x(t) as h — > with nh = t). The order of a convergent 
/c-step method cannot be higher than k + 2. 

The orbital instability of the Stdrmer-Cowell methods can be explained by a simple example 
(Lambert and Watson 1976). Consider equation (|j), whose general solution is the periodic function 
x(t) = Acos(ujt) + Bsin(u>t). Applying the multistep method (||) to this equation we obtain the 
difference equation 

k 

+ H 2 (3i)x n+i = 0, H = hu, (10) 

whose general solution is 

k 

x n = J2AjZ?, (11) 

3=1 

where the Zj (j =1, 2, . . . , k) are the roots, assumed distinct, of the polynomial 

Q(Z;H 2 ) = p{Z) +H 2 a{Z). (12) 

The roots Zj of O are perturbations of the roots Zj of p\ let Z\ and Z2 be the perturbations of the 
principal roots. The problem with the Stormer-Cowell methods is that, if k > 2, the roots Z\ and 
Z2 do not lie on the unit circle. No matter how small H 2 is the orbit grows or shrinks, and the 
energy error increases linearly with time. 

Lambert and Watson (1976) showed that this orbital instability can be avoided if the multistep 
coefficients are chosen to be symmetric, as in equation (||). A multistep method is said to have an 
interval of periodicity (0, Hq ) if for all < H < Hq the roots of Q(Z; H 2 ) satisfy 

|Zi| = |Z 2 | = 1, \Zj\<l (j = 3,...,k). (13) 
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Lambert and Watson proved that a convergent multistep method with a non-zero interval of peri- 
odicity must be a symmetric method and must have an even order. For a symmetric method the < 
in can be replaced by =, and hence inside the interval of periodicity the roots all lie on the unit 
circle. The solution (|ll]) of equation (|j) is then guaranteed to be periodic (or quasi-periodic). 

The requirement that a /c-step symmetric multistep method have order k (for an explicit method) 
or k + 2 (for an implicit method) does not determine the a and (3 coefficients uniquely when k > 2, 
because for a symmetric method the equations Cj = with j odd are not independent of the 
equations with j even. There are thus some free coefficients that can be chosen by the user, although 
their range is restricted by the requirement of zero-stability. Lambert and Watson (1976) gave 
examples of explicit methods with orders 2, 4, and 6, and implicit methods with orders 4, 6, and 
8. Quinlan and Tremaine (1990) gave examples of explicit methods with orders 8, 10, 12, and 
14. Table 1 lists the coefficients of five explicit methods that will be discussed in what follows. 
The methods SY8, SY10, and SY12 are the 8th-, 10th-, and 12th-order methods of Quinlan and 
Tremaine; SY8A and SY8B are 8th-order methods that have not previously been published. Table 2 
lists the spurious roots of p{z) for these five methods. 
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Table 1: Coefficients of the symmetric multistep methods (only half are listed; the others are given by 
OLi = a k -i, Pi = Pk-i)- 
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6.000 


60.00 


SY8A 
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8.000 


16.00 
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4.500, 6.000, 9.000 
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Table 2: Spurious Roots of p(z) for the methods in Table [l| 



3 Origin of the resonances and instabilities 

The origin of the resonances and instabilities will be explained using a Kepler orbit as an example. 
The resonances and instabilities occur for all nonlinear oscillatory problems, not just for the Kepler 
problem; other examples will be given later. 
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3.1 A simple explanation for the instabilities 

We start with a simple explanation for the instability that occurs when a circular orbit is integrated 
with the method SY8 using 60 steps per orbit. The spurious roots of p(z) for the method SY8 are 
located on the unit circle at angles ±4vr/5, ±2tt/5, and ±2vr/6 (i.e., at ±144°, ±72°, and ±60°). At 
60 steps per orbit the spurious roots of Q(Z;to 2 h 2 ) differ little from those of p(z); the difference will 
be ignored here. It is the roots at 2tt/5 and 2-7r/6 (or —2tt/5 and — 2tt/6) that cause the trouble. 
The root at 2ir/5 allows a 5-step oscillation: in the absence of any forces, equation (||) admits a 
solution x n = exp(27rin/5). Similarly, the root at 27r/6 allows a 6-step oscillation. By themselves 
the oscillations would be harmless for this problem as long as the start-up routine is accurate. The 
trouble arises when the 5- and 6-step oscillations can resonate. 

Consider first the 5-step oscillation. The perturbations to the x and y coordinates can both 
oscillate with a 5-step period. If the y oscillation is 90° ahead of the x oscillation, the perturbation 
goes around in a counter-clockwise sense with a 5-step period. Assume that the orbit also moves in 
a counter-clockwise sense. During one orbital period the perturbation goes around 60/5 = 12 times. 
Because the perturbation goes around in the same sense as the orbit, however, the perturbation 
to the orbital radius completes only 12 — 1 = 11 oscillations. Now consider the 6-step oscillation. 
This time assume that the y oscillation is 90° behind the x oscillation, so that the perturbation goes 
around in a clockwise sense with a period of 6 steps. During one orbital period the perturbation 
goes around 60/6 = 10 times, but the radial perturbation completes 10 + 1 = 11 oscillations. This 
leads to the resonance: the radial perturbations from both the 5-step and 6-step oscillations can 
go around 11 times in one orbital period. The perturbation analysis given below shows that the 
resonance causes an instability. 

The explanation was verified by checking that the instability can be enhanced by adding noise 
to the initial conditions with the right frequency and phase. When noise was added with a 5-step 
component polarized in a counter-clockwise sense and a 6-step component polarized in clockwise 
sense, the instability became obvious sooner, but when the polarizations were reversed, so that the 
5-step component was clockwise and the 6-step component counter-clockwise, the instability was 
delayed. 

The explanation predicts that instabilities will occur for other symmetric multistep methods at 
stepsizes at which the number of steps per orbit N satisfies 

N +1 = ^-1, (14) 



2tt/0/ 2ir/0j 

where exp(i^) and exp(i0j) are spurious roots of p(z). The order of the method must be at least six 
for an instability of this type to occur, since the polynomial p(z) must have at least four spurious 
roots. According to the prediction the method SY10 of Quinlan and Tremaine should be unstable 
for circular orbits at 84 steps per orbit, and the method SY12 should be stable as long as the number 
of steps per orbit is at least 36; both predictions were verified, along with similar predictions for 
other symmetric multistep methods. 



For eccentric Kepler orbits equation (14) must be generalized to allow for the different frequency 
components present in the motion. The forces in an eccentric Kepler orbit (with unit major axis) 
can be expanded as (Kovalevsky 1967) 



r(t) 3 

y(t) 

r(t) 3 



9 [ J q+i(Q e ) ~ Jq-i(qe)] cos{qut), (15) 



9=1 



Vl - e 2 Y Q [Jq+i(Qe) + J q -i(qe)] sin(gwt), (16) 

9=1 
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where the J q are Bessel functions of the first kind. The forces can be built up from components with 
frequencies that are integral multiples of the fundamental frequency u. The "1" on the left-hand side 



of equation (14) must therefore be allowed to take on the integer values 1, 2, 3, . . . , corresponding to 
the fundamental frequency and the higher harmonics, and similarly the "1" on the right-hand side 
must be allowed to take on the same values, independently of the value assumed on the left-hand 
side. This leads to the prediction of instabilities at N = 60, 90, 120, 150, etc. The width of the 
instability at small eccentricities (when plotted as in Figure is found to vary with the eccentricity 
as ~ e° at N = 60, as ~ e 1 at N = 90, as ~ e 2 at N = 120, and so on, as expected since the different 
frequency components have amplitudes that vary as powers of the eccentricity. 

3.2 Perturbation analysis for a circular orbit 

A perturbation analysis can be used to show that circular orbits are unstable when the condition ( |i~4| ) 
is satisfied. Consider a circular orbit of radius R in a two-dimensional axisymmetric potential (j>{r). 
The circular frequency u is 

uj( r ) 2 = 4>'{r)/r. (17) 

Let x n and y n be the x and y coordinates at the nth time step. Define z = x + iy, r = \z\, and write 
equation @ as 

k k 

zZ a J Z n+j = h2 Yl Pj F n+j, (18) 
3=0 3=0 

where the complex force F is 



F = U + ify 



-cf>>(r). 



(19) 



Provided that the stepsize h lies within the interval of periodicity, and that the radius is assumed 
to be fixed, so that equation (^) is linear in z, the principal roots Z^ 1 (in fact all the roots) of (|i~8| ) 
will lie on the unit circle, and the numerical solution will have the form 



rz: 



Z p ~e 



iuih 



(20) 



where Z p is the principal root of Q(Z; u) 2 h 2 ) corresponding to the assumed counter-clockwise rotation. 
Now consider a perturbed circular orbit 



z n = RZ£(l+u n ), 

where RZJ}u n is a small perturbation at time step n. The perturbed radius is 



(21) 



Rn 



(ZnZ*) 1 / 2 ~i?[l + K + <)/2], 



(22) 



where the * denotes complex conjugation. Substituting (21) into equation ( |lg| ) and linearizing the 
resulting equation we find 



3=0 



where 



UJ-2 



3=0 



-<f>'(r) + <t>"{r) 

T 

-0'(r) - <f>"(r) 
r 



(k 2 - 2uj' 



(23) 

(24) 
(25) 
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and where k is the epicyclic frequency, given by (Binney and Tremaine 1987) 

«2( r ) = r^L + Au 2 = cp"(r) + V(r). (26) 
ar r 

A trial perturbation of the form 

u n = A l S n + A 2 S* n (27) 

leads to a solution for 5 that is independent of the step number n if the following two equations are 
satisfied: 

k k 

A^fajZlS* + co 1 h 2 f3 J Z^) = AZ^untfPjZiS*; (28) 

3=0 3=0 
k k 

A 2 J2{a3 Z iS* j +u l h 2 /3 j ZiS* j ) = A\Y,"2 h2 Pi z i s * j - ( 29 ) 

3=0 j=0 

Taking the complex conjugate of the second equation, we get two equations for the two unknowns 
A\ and A\. The determinant of the system must vanish, which gives (using H = uoh) 

D(S) = fl(SZ p ; H 2 )Q(S/Z P ; H 2 ) - uj 2 h 2 Wl(SZ p ; H 2 )a(S/Z p ) + Q(S/Z p ; H 2 )a(SZ p )} = 0. (30) 

D(S) is a symmetric polynomial of degree 2k in 5, with real coefficients. There are several things 
to note about the roots of this polynomial. If 5 is a root then so are 5* , 1/5, and 1/5*. 5 = 1 
is always a root, since Q(Z p ;H 2 ) = fl(l/Z p ; H 2 ) = 0. The root 5 = 1 corresponds to unperturbed 
motion, because for this root A\ = —A\ and hence the linearized perturbation to the radius is zero. 
If h = then D(S) = p(S) 2 , whose roots are the same as those of p(z), with each root Si having 
twice the multiplicity of the corresponding root Zi of p(z). Thus when h = the spurious roots of 
p(z) (assumed to be distinct) are double roots of D(S), and 5 = 1 is a root of multiplicity four. If h is 
small but nonzero the roots that were double roots of D(S) when h = split, and are approximately 
ZjZi^ 1 (j = 3, . . . , k), where the Zj are the spurious roots of Q(Z; H 2 ). The root 5 = 1 is a double 

root when h ^ 0, and two new roots appear at approximately Z p K ^ ', corresponding to the usual 
epicyclic oscillations with frequency u> ± n. 

As h increases away from the roots of D(S) move around the unit circle as just described. 
Problems arise, however, near stepsizes where Z 2 Z\/Zj = 1 for some spurious roots Z\ and Zj. 
When this happens the approximate solutions for two of the roots coincide, Z\Z p = Zj/Z p , and a 
more careful analysis is needed. The troublesome stepsizes can be estimated by replacing Z p by 
ex.p(iujh) and Zi and Zj by the corresponding roots of p(z), which we write as z\ = exp(i#/) and 
Zj = exp(i9j). The troublesome stepsizes then occur when 

6j -6i = 2uh = 4tt/N, (31) 

where ./V is the number of steps per orbit. This is the instability criterion ( |l4] ) that was given earlier. 

Figure |3| shows the maximum magnitude of the roots of D(S) for the method SY8 used near 60 
steps per orbit with a Kepler potential <f>(r) = — 1/r. The solid line is the exact solution found by 
numerical calculation of the roots. The dashed line is an approximate solution found by simplifying 
D(S) by writing 

8 

UiSZ} 1 ; H 2 ) ~ [] (SZp 1 ~ Z m , r ), (32) 

m=l 

where the Z m)T are the roots of £1 at the resonant stepsize (N ~ 60.455), and by replacing SZ p 
by Zi tr and S/Z p by Z^ r everywhere except in the terms (SZ p — Z^ r ) and (S/Z p — Zj ir ); the same 
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Figure 3: Maximum value of |Sj| (i = 1, . . . , 16) for the roots of D(S) using the method SY8 with a Kepler 
potential <j>{r) = — 1/r. The solid line is the exact result from a numerical calculation of the roots; the 
long-dashed line is the approximate solution described in the text. The short-dashed line marks the value 
TV ~ 60.455 at which Z^Zi/Zj = 1 for two of the spurious roots of 0. 

replacements are made in the a functions. We thus obtain from D(S) = a quadratic equation 
in S, whose roots are easily found. The exact and approximate solutions both show that there are 
roots of D(S) that lie off the unit circle when N is in the range 59.2-60.4. The amount by which the 
roots move off the unit circle depends on the value of k 2 — Auj 2 . For a harmonic oscillator potential 
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<f>(r) ~ r 2 the roots of D(S) do not move off the unit circle, since k 2 — 4u; 2 =0, and the instability 
does not occur. 



3.3 Location of the resonances 

The resonances in Figure [I] are easier to predict than the instabilities. Since the force acting on an 
eccentric Kepler orbit is a superposition of components with frequencies that are integral multiples 
of the fundamental frequency uj, the resonances occur when the principal root Z p ~ exp(iujh) raised 
to some integral power q coincides with one of the spurious roots Zj, i.e., when 

N~^-, g = l,2,3,..., (33) 

where Zj = exp(i6j) is a spurious root of p(z). This correctly predicts the resonances for the 
method SY8 at stepnumbers that are multiples of 5 and 6 (and also 2.5, although these resonances 
are usually too weak to be seen). The amplitude and growth rate of the resonance decrease as q 
increases, because for small eccentricities the Bessel functions in equations (|l5|) and ((l^) decrease 
rapidly with q. 



4 Further examples 

Three more examples will be given to show how the resonance and instability locations depend on 
the multistep method and on the equation being integrated. 



4.1 Kepler orbits with a 12th-order symmetric method 

The 12th-order symmetric method SY12 is expected to be better than the 8th-order method SY8 
for circular Kepler orbits, since the 12th-order method is free of instabilities for these orbits as 
long as the integration takes at least 36 steps per orbit. This is confirmed in panel (a) of Figure || 
(the energy error of the symmetric method is caused by roundoff error at most stepsizes in this 
panel). For eccentric orbits, however, the 12th-order method can be as troublesome as the 8th-order 
method, if not more so, because the extra spurious roots allow more opportunities for resonances and 
instabilities to occur, and because the spurious root at 2ir/9 on the unit circle leads to a resonance 
that can be excited even at low eccentricities. Panel (b) of Figure || shows resonances at stepnumbers 
that are multiples of 4.5, 6, and 9, and instabilities at stepnumbers that are multiples of 36, although 
these are not the only instabilities; there are others in the figure, and more would be present if the 
stepsizes had been sampled more finely. 



4.2 Orbits in a logarithmic potential 

The Kepler potential is special in that all bound orbits are closed and have uj = k, i.e., the azimuthal 
period T a = lis juj is equal to the radial period T r = 2tt/k (equation ([26|) gives k for a circular orbit). 
For most potentials T a ^ T r , and the analysis given previously must be modified. Consider an 
eccentric orbit in a logarithmic potential 4>{r) = log(r), for which (in the limit of a circular orbit) 
k = y/2u). A Fourier analysis of the motion contains frequencies uj + qn (where q is an integer), and 
hence the prediction (|33"|) for the resonant stepsizes must be modified to read 



N ^( l + ^) = ^( 1+ -^.), , = 0,1.2,..., (34) 



uj J 9j \ T r 
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Figure 4: Maximum fractional energy error during integrations of a Kepler orbit for 25000 periods using 
the method SY12 (solid lines) and the 12th-order Stormer method (dashed lines) with about 1600 different 
stepsizes, plotted as a function of the number of steps per orbit: (a) circular orbit; (b) eccentric orbit (e = 0.2). 

where N = T a /h is the number of steps per azimuthal period. This prediction is verified for the 
method SY8 in Figure ||, which shows results from integrating an orbit with initial conditions xq = 1, 
yo = 0, x'o = 0, y'o = 1.1, for which T a /T r = 1.41536, close to the value v2 for a circular orbit (the 
stepsizes were not sampled as finely in this figure as in the other figures). There is an instability 
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Figure 5: Maximum fractional energy error during integrations of an eccentric orbit in a logarithmic potential 
for 10000 periods with the method SY8, using 500 equally-spaced stepsizes (see text). 

at 60 steps per azimuthal period, just as with the Kepler potential, and resonances at the N values 
predicted by equation (|34|), taking 9j to be 2ir/5 (the short dashed lines in Figure ||) or 2n/6 (the 
long dashed lines). 
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4.3 The one-dimensional hard spring 

As a final example we consider the one-dimensional hard-spring equation 

ar"(i) = - x (t) 3 . (35) 

While this is not an orbital equation, the motion is similar to an eccentric orbit in that it is periodic 
and contains frequencies that are integral multiples of the fundamental frequency, although in this 
case only the odd multiples are present. Figure [6| shows the errors obtained by integrating this 
equation with the method SY8. As expected, there are resonances at stepnumbers that are odd 
multiples of 5 or 6, such as 65, 66, 75, 78, etc.; the even multiples are missing because of the missing 
frequencies in the motion. There are instabilities at 60 steps per orbit and at integer multiples of this 
number: 120, 180, etc.; the instabilities at 90, 150, etc., are missing, again because of the missing 
frequencies. 

5 A search for better symmetric methods 

A question that naturally arises is whether Quinlan and Tremaine (1990) could have chosen better 
multistep coefficients to reduce the resonance and instability problems. The answer is yes for the 
8th- and lOth-order methods, but probably no for the 12th-order method. 

There are three properties that we would like a symmetric multistep method to have: 

1. The method should have a large interval of periodicity. 

2. To avoid instabilities the spurious roots of p{z) should be well spread out on the unit circle. 

3. To avoid resonances the spurious roots of p{z) should be far from z = 1. 

Fukushima (1998) has searched for multistep methods with large intervals of periodicity. It is 
properties |2| and |3| that concern us here. Unfortunately these properties are not compatible: the 
spurious roots cannot be well spread out on the unit circle and at the same time be far from z = 1. 
A compromise must be made, which becomes more difficult the more spurious roots there are, i.e., 
the higher the order of the method. 

A systematic search was made through symmetric multistep methods with a coefficients drawn 
from the set {0, ±1/8, ±2/8,. . . ,±7/8,±l,±2,±4} (with at = 1). This set was chosen because the 
integration method suffers less from roundoff error if the a's are integral powers of 2; it is unlikely 
that much better methods would have been found by letting the a's take on arbitrary values. For 
each set of coefficients the roots of p(z) were computed and checked to see if and where they lie 
on the unit circle. The most promising methods (based on the three properties given above) were 
compared in integrations of eccentric Kepler orbits. 

Consider first the 8th-order methods. For any even integration order k, the choice 1, —2, +2, —2, 
+2, . . . , for the a coefficients spreads the spurious roots on the unit circle as evenly as possible, and 
tends to minimize the instabilities. That choice was made for the method SY8A listed in Table |], 
which has the largest interval of periodicity of the 8th-order methods that were tested. This method 
is much more stable than the method SY8; the integration results in Figure ^ show no signs of 
instabilities even for an orbit with an eccentricity e = 0.3. The drawback of the method SY8A, 
however, is the spurious root at 2ir/8, or 45°, which leads to resonances (at stepnumbers N that are 
multiples of 8) that are stronger than for SY8. The method SY8A is certainly an improvement over 
SY8, but the resonances are annoying. 

The resonances can be reduced by picking a method whose spurious roots are farther from z = 1. 
An example is the method SY8B, whose spurious root closest to z = 1 is at 76.96°. (It is easy 
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Figure 6: Maximum energy error measured during integrations of the hard-spring equation (35) for 50000 
periods starting from the initial conditions xo=0 and x' = 1/V2, using the method SY8 with 1600 different 
stepsizes. 



to find methods whose resonances are weaker than those of SY8B, but these methods have much 
poorer stability properties.) Figure shows that the resonances are much weaker with this method 
than with SY8A and SY8: they fall off faster as the number of steps per orbit is increased, being 
almost unnoticeable when the e = 0.15 orbit is integrated with more than 75 steps per orbit, and the 
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Figure 7: Error versus stepsize during integrations of Kepler orbits for 25000 periods with the 8th-order 
symmetric methods SY8A and SY8B, using 1500 different stepsizes. 

e = 0.30 orbit with more than 120 steps per orbit. But the price we pay for this gain is a decrease 
in stability, as revealed by the spikes in the error plot, especially at the higher eccentricity, e = 0.3. 
The instabilities can be avoided if the stepsize is chosen small enough, however, and are not as bad 
as with the method SY8. The method SY8B is the most promising of the 8th-order methods that 
were tested. 
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The search for better 10th- and 12th-order symmetric methods was not as successful. The method 
SY10 of Quinlan and Tremaine (1990) is unstable at 84 steps per orbit, even for a circular orbit. 
The method SY10A (like SY8A, with the spurious roots spread out evenly on the unit circle) is 
much more stable, but has an annoying resonance at stepnumbers N that are multiples of 10. None 
of the other lOth-order methods tested was much better than SY10A. For the 12th-order methods 
none was found to be better than SY12. The method SY12A (again, with the spurious roots spread 
out evenly on the unit circle) has an annoying resonance at stepnumbers that are multiples of 12, 
and its stability properties are no better than those of SY12. 

Fukushima (1999) has tested some implicit symmetric methods to see if they are affected by the 
resonance and instability problems as badly as the explicit methods. The disadvantage of implicit 
symmetric methods is that the corrector step must be applied at least twice for the benefits of the 
symmetry to be realized, so that for the same stepsize an implicit method requires at least three 
times the number of force evaluations required by an explicit method. Unless implicit methods can 
be found with much weaker resonances and instabilities than the explicit methods, this extra work 
would be better spent using an explicit method with a smaller stepsize. 

6 The use of symmetric methods for integrating planetary orbits 

The examples considered so far have been chosen for their simplicity and pedagogical value. We 
now test the method SY12 on a more complex example, typical of those encountered in research 
problems, to get a better idea of how the method compares with a high-order Stormer method and 
how its effectiveness is reduced by the resonances and instabilities. The example is the long-term 
integration of a planetary system, the problem that was the motivation for the high-order symmetric 
multistep methods of Quinlan and Tremaine (1990), and a research problem on which the method 
SY12 has been used (the 3 Myr integration of all nine planets by Quinn, Tremaine, and Duncan 
1991). We consider an idealized solar system containing only two planets, Jupiter and Saturn, as 
this is sufficiently complex to suggest how the method will work for a system with more than two 
planets. 

The initial conditions for Jupiter and Saturn are taken from Standish (1990). The initial major 
axes are aj ~ 5.203 AU and as — 19.280 AU, giving orbital periods Pj ~ 4332.8 and Ps — 30905 (in 
days, the unit of time in this discussion). The initial eccentricities are ej ~ 0.048 and es — 0.051. 
The orbits were integrated for 1 Myr using 1000 different stepsizes h spaced equally in 1/h between 
h = 50 and h = 81, corresponding to approximately 86.6 and 53.5 steps per orbit for Jupiter. 
To speed the calculations the energy error was measured every 5th integration step, rather than 
every step. Jupiter's longitude at the end of the integration was compared with an accurate value 
determined by an integration with a much smaller stepsize {h = 10). The integration errors are 
plotted versus stepsize in Figure ||. 

The two planets were first integrated separately, with no gravitational interaction between them. 
The maximum energy errors are shown in panel (a). The errors for Saturn in this panel are caused 
by roundoff error and are not interesting. The Stormer method is unstable if it is used for Jupiter's 
orbit with a stepsize larger than h ~ 57. The symmetric method is stable for Jupiter's orbit with 
stepsizes as large as h ~ 80, although the effects of resonances are seen when the number of steps 
per orbit is a multiple of 9 (near the h values 53.5, 60.2, and 68.8, corresponding to 81, 72, and 63 
steps per orbit), and also, to a lesser extent, a multiple of 6 (near h = 72.2, or 60 steps per orbit). 
The spike in the error near h = 80 (54 steps per orbit) is an instability, not a resonance, as there 
is a sudden growth in the error to a large value early in the integration. Away from the resonances 
the maximum energy error is more than 100 times smaller with the symmetric method than with 
the Stormer method. 
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Figure 8: Error versus stepsize for 1 Myr integrations of the Jupiter-Saturn system using the 12th-order 
symmetric method SY12 (solid lines) and 13th-order Stormer method (dashed lines). The ordinate gives 
the common logarithm of the error. Panel (a) shows the maximum energy error in the orbits of Jupiter 
(upper lines) and Saturn (lower lines) when the planets are integrated separately, with no interaction between 
them. The other three panels include the gravitational attraction between Jupiter and Saturn, and show the 
maximum (b) and average (c) energy errors during the integration, and (d) the longitude error for Jupiter at 
the end of the integration. 
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The planets were then integrated together including their gravitational interaction, in the hope 
that it might detune the resonances and remove the spikes from the error plot. This is not what 
happened. Panel (b) shows a profusion of new resonances resulting from the high-frequency pertur- 
bation terms in the Jupiter-Saturn interaction. There is still the instability near h = 80 that was 
present for Jupiter alone, plus a narrow instability near h = 78.55 that was not present for Jupiter 
alone, but it is hard to identify the resonances at 60, 63, 72, and 81 steps per Jupiter orbit, as they 
are lost in the other resonances. 

With the symmetric method the maximum energy error away from the resonance peaks is no- 
ticeably larger for the Jupiter-Saturn system than for the single-planet (Jupiter) system, whereas 
for the Stormer method the errors for the two systems are comparable. At first this suggests that 
the symmetric method is not much better than the Stormer method for the two planet system. But 
panels (c) and (d) suggest the opposite conclusion. With the symmetric method the average energy 
error is several orders of magnitude smaller than the maximum error, whereas with the Stormer 
method the average and maximum errors are comparable. The error in Jupiter's longitude at the 
end of the integration is much smaller with the symmetric method than with the Stormer method. 
Even at the resonance peaks the longitude error is more than an order of magnitude smaller with 
the symmetric method than with the Stormer method, and away from the resonances the difference 
is more than three orders of magnitude. The narrowness of the resonances is difficult to appreciate 
from the figure. Out of the 310 stepsizes that were tested in the range h = 60-70, only 16 (or 8) 
resulted in longitude errors larger than 10~ 5 (or 10~ 4 ). 

These results show that the maximum energy error can be an unduly strict criterion to use for 
evaluating the performance of a symmetric method. In a planetary integration the position errors 
are determined by the average energy error, which for the symmetric methods is much smaller than 
the maximum error. The maximum error is important for calculations of instantaneous orbital 
elements that depend on the velocity, such as the major axis, eccentricity, and inclination. In the 
3 Myr integration of the solar system by Quinn et al. (1991) these elements were digitally filtered 
to remove oscillations with periods smaller than 500 yr; the filtering probably removed any spurious 
oscillations caused by the symmetric method, and the errors in the output elements were probably 
closer to the average errors than to the maximum errors. The comparison that Quinn et al. made 
of their results with the results of a shorter but more accurate Stormer integration (with a smaller 
stepsize) showed satisfactory agreement. 

Despite the resonances and instabilities, then, symmetric methods can still be a better choice 
than Stormer methods for long integrations of planetary orbits provided that the user is aware of 
the dangers. 
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